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ABSTRACT 


'rtiis  report  concerns  the  prediction  of  the  flow 
characteristics  of  an  isothermal  free  jet.  A computer  program 
has  been  developed  similar  to  that  of  Spalding  and  Patankar. 
Reynolds  stress  equations  are  used  so  that  not  only  turbulent 
shearing  stress,  but  also  turbulent  kinetic  energy  and 
dissipation  can  be  calculated.  This  program  is  rather  short, 
about  280  statements,  and  for  a moderate  number  of  points 
(usually  about  15) , requires  only  five  seconds  per  run  for  the 
Amdahl  47P/V6  computer.  The  results  coqpare  fairly  well  with 
experiments  in  two-dimensional  as  well  as  in  axi-symmetric  jets. 
It  is  found  that  the  similarity  assuqption  is  only  approximate. 
Also  the  results  can  be  somewhat  different  for  different  initial 
input  turbulence  conditions.  Therefore,  to  compare  the 
experimental  results  and  to  interpret  their  accuracy, 
particularly  when  no  detailed  measurements  are  made  at  the  jet 
orifice,  should  be  done  cautiously.  The  variations  due  to  the 
assigned  constants  in  the  closure  model  are  also  briefly 
discussed. 
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INTRODUCTION 


Combustion  in  dump  combustors  is  a complex  process  which 
involves  mixing,  mass  transfer,  chemical  reaction  as  well  as 
circulations.  The  investigation  using  the  Reynolds  stress 
equation  for  predicting  the  flow  characteristics  of  an  isothermal 
free  jet  is  one  of  the  very  first  steps  toward  understanding  the 
mixing  process.  Based  on  the  boundary  layer  approximation  of  the 
jet  mixing  and  neglecting  the  effects  of  a wall,  the  coupled 
equations  involving  momentum,  turbulent  shear  stress,  kinetic 
energy  and  dissipation  are  much  simplified  with  a proper  closure 
model.  The  numerical  prediction  by  the  Spalding  method  with  the 
von  Mises  transformation  is  investigated  with  some  modification 
for  the  free  jet  calculation.  A computer  program  is  developed  for 
use  such  that  it  is  relatively  easy  for  the  user  to  read  and  to 
modify  the  program  for  his  needs.  Further  development  should 
include  the  chemical  species  such  that  free  jet  combustion  can  be 
predicted. 

Since  the  APOSR-IFP-STANFORD  Conference  1968  the  prediction 
of  the  turbulent  flow  has  shifted  its  emphasis  to  turbulent  field 
methods  which  involve  the  Reynolds  shear  stress  turbulent  energy 
or  turbulent  dissipation  with  varying  closure  models  as  discussed 
in  that  report  [1].  In  1970,  Reynolds  [2]  presented  a brief 
survey  of  the  state  of  the  art  for  computation  of  turbulent 
flows.  Herring  and  Mellor  ( 3 J developed  a computer  program  which 


is  used  fairly  widely  in  the  United  States.  At  Imperial  College, 
Spalding  and  his  coworker  [4]  worked  over  ten  years  on  a program 
which  is  large  and  versatile  with  different  turbulence  models. 

His  program  is  used  in  Europe  as  well  as  in  the  U.S.  Because 
the  program  is  large  and  versatile  it  requires  some  training 
and  fluid  mechanics  background  to  use  the  program  properly. 

If  the  user  wants  to  modify  the  program  for  his  needs,  he  usually 
encounters  many  difficulties.  It  is  well  known  that  it  is 
difficult  to  read  and  modify  a computer  program,  especially  if  it 
is  a large  and  complicated  one,  unless  the  user  is  quite  familiar 
with  the  detailed  procedure.  Launder  [S,6]  used  the  Reynolds 
stress  equation  in  a sequence  of  his  papers  on  predicting 
turbulent  flow.  His  criticism  is  the  common  one  i.e.  there  are 
too  many  constants  and  it  is,  in  a sense,  a conplicated  curve 
fitting  technique.  Nevertheless,  Launder 's  approach  does  give 
fairly  good  results  in  general.  However,  there  are  no  reports  on 
the  predictions  of  axi-symmetr ic  jets  which  involve  turbulent 
shear  stress  and  kinetic  energy  as  well  as  dissipation.  This 
report  is  for  such  an  investigation. 
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GOVERNING  EQUATIONS 


1.  Turbulence  Closure  Model 

For  a fluid  of  uniform  density  p the  set  of  differential 
equations  governing  the  transport  process  can  be  written  in  the 
following  form: 

Equation  of  continuity 


9U, 

dx; 
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Equations  of  momentum 


Equations  for  (kinematic)  Reynolds  stress 
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Equation  for  dissipation 


= /-*££**  £££*k 

ZD*  axA  ( 9^  9ZX  9*  *x4  ) 9XX  9*4 


-JiJL  92^  > 


where  U is  the  mean  flow  velocity  of  the  main  flow, 


ux»  is  the  fluctuation  velocity,  i.e.  u^  * 0, 


^ , d2t;  dll; 

£ — ~~2/  — is  the  dissipation  and 


£ is  the  dissipation  fluctuation 

Since  there  are  more  unknowns  than  the  number  of  equations, 
some  closure  assumptions  have  to  be  made  in  order  to  solve  the 
equations.  Based  on  the  information  on  isotropic  turbulence, 
homogenous  turbulence  and  pure  shear  flow  for  large  Reynolds 
nunber,  Hanjalic  and  Launder  [5]  propose  the  following  form  of 
closure  for  Reynolds  stress  and  dissipation. 
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where  ■ 2.8  and  C2  ■ 0.070^ 

There  are  six  constants  involved  in  the  closure  model  as 
tabulated.  Thus  we  have  five  equations  with  five  unknowns,  0,  V, 
uv,  k and  £ • 
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This  set  of  equations  is  parabolic  and  the  initial 
conditions  for  U,  uv,  k and  £ are  usually  not  available  at  the 
beginning  station,  hence,  some  guess  work  of  their  distributions 
must  be  made.  It  is  observed  that  owing  to  the  uncertainty  of  £ 
and  k in  the  earlier  stages,  the  range  of  the  ratio  k/e  might  be 
quite  large  by  simply  assuming  some  arbitrary  distributions  of 
and  k respectively.  One  way  to  overcome  this  is  by  assuming  that 
the  kinetic  energy  k is  proportional  to  the  dissipation  £ across 
the  stream.  This  idea  comes  from  the  limiting  case  that  as  £ 
approaches  zero  or  k approaches  zero,  the  ratio  k/£  trust  approach 
a finite  value;  otherwise  an  artificial  singularity  is 
introduced.  This  assumption  seems  particularly  appropriate  for 
free  boundary  flow.  From  dimensional  considerations,  it  is  found 
that  k/£  is  proportional  to  Qs/t/  for  free  jet  flow  where  SoS  is 
the  conventional  half  width  where  the  velocity  is  half  of  the  center 
velocity  Uc  . Thus,  it  appears  to  be  proper  to  introduce 

fa  = £ =cKCc + co 


which  gives  the  asynptotic  expression 


^ DC  3/X 


for  the  plane  jet  and 


for  an  axi-synmetric  jet. 
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Ttie  constant  C represents  the  adjustment  of  the  virtual 
origin.  If  this  approach  is  adopted,  it  serves  the  following 
advantages  in  the  final  form  of  the  equations. 

(i)  The  equation  for  dissipation  is  eliminated  at  the 
beginnning  stages,  thus  reducing  the  system  of  five  differential 
equations  to  four 

(ii)  Only  initial  profiles  of  U,  uv  and  k are  needed. 

(iii)  The  sensitive  constant  C_  , is  avoided  (at  least  at 

/ 

the  beginning  stages).  As  reported  by  Launder  (6),  in  the 
dissipation  equation  is  such  a sensitive  constant  that  a small 
change  of  C*  , will  result  in  a large  change  in  the  predictions. 

As  a matter  of  fact,  launder  later  proposed  to  use  the  value  1.44 
as  compared  to  1.4S  as  had  been  suggested  previously. 

(iv)  It  also  provides  the  extension  to  calculate  the 
dissipation  from  the  original  five  equations  by  using 

f(x)  - C (C+  8 /u)  for  the  development  of  the  initial  stages  and 

JC  o%€*  c 

then  use  f(x,y)  * k/£  directly  for  further  downstream  calculation. 
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NUMERICAL  METHOD  AND  CALCULATION  PROCEDURE 

1.  Transformation  from  physical  coordinates  (xfy)  to  streamline 
coordinates  (x,^) 

The  afore-mentioned  five  equations  for  U,  V,  uv,  k and  <5 

can  be  reduced  further  to  four  equations  for  U,  uv,  k and  £ by 

the  von  Mises  transformation  from  physical  coordinates  (x,y)  to 

streamline  coordinates  (x,  'iff ) as  discussed  by  Pantanka  & 

Spalding  [4] . Introducing  the  stream  function  'i/T  and  a 

T fr 

nondimensionalized  stream  function  uJsrr?  where  "m  represents  the 

2tf 

streamline  at  the  edge  of  the  flow.  Thus 

Finally  it  can  be  arranged  in  the  following  form  in  0 
coordinates. 
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It  is  noted  that  the  equations  are  linearized  so  that  the  matrix 
of  the  coefficients  is  in  tri-diagonal  form.  The  solution  of 
the  equations  can  be  put  in  the  following  form  as  discussed  in 
Schlichting  (7)  by  letting 

= * a-  m 

J \J  "ZV*/  J 

where 


a.  = % onj  h.  = £ — us> 

The  boundary  condition  for  free  jet  by  using  symmetry  at  centerline 


y ■ 0: 

rj,/  J 

■ — = 0 gives 

ay 

G / " 0 

H,  = 1 

uv  * 0 

gives  G(  * 0 

Hy  = 0 

&A  „ n 

gives  G,  * 0 

»,  * P 

and  "for- 

3£  / 

t=  0 o/sO 

ay 

Thus  after  calculation  of  G • and  H,*  from  the  center  toward  the 

J v 

boundary  edge  for  j * 2,3....  N+l  the  unknowns  <3?,  are  successively 

"V 

obtained  from  the  edge  toward  the  center  by  Equation  (121. 

3.  Entrainment 


In  the  marching  forward  calculation,  the  non-dimensional 
stream  function  (J  is  constant  and  at  the  edge  CO  * 1.  However, 
the  stream  function  at  the  edge  is  a function  of  x,  and  it  is 
expected  that  will  increase  as  x increases  due  to  entrainment 
By  arbitrarily  choosing  a point  close  to  the  edge  of  the 

ctx. 

boundary  and  arbitrarily  requiring  a value  of  the  velocity  at  CO^  , 
say  Ug  ■ fraction  of  the  center  velocity  Uf,  the  entrainment  m 
can  be  evaluated  from  Equation  (7).  Alternatively,  the  entrainment 
rate  is  properly  chosen  as  is  the  case  in  this  program. 
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The  computer  program  has  been  written  for  a free  turbulent  jet 

with  output  in  decay  of  center  velocity,  velocity  profile, 

turbulent  shear  stress,  turbulent  kinetic  energy  and  dissipation. 

The  goal  was  to  write  the  program  concisely  so  that  it  is 

easier  for  the  user  to  read  and  modify  the  program  for  his 

needs.  This  program  has  been  tested  for  sensitivity  in  initial 

conditions,  constants  and  entrainment  rate.  In  the  program,  the 

total  momentum  (FLUX)  is  calculated  and  is  compared  with  initial 

total  momentum  (SM) , the  difference  DM=  FLUX  - SM  serves  a check 

for  conservation  of  momentum.  Usually  the  ratio  is  about  a 

SM 

few  percent. 
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RESULTS  AND  DISCUSSIONS 

Based  on  the  cepoct  by  Tsuei  [8),  which  has  been  tested  on 
laminar  and  turbulent  two  dimensional  jets,  this  computer  program 
is  an  extension  to  the  case  of  the  axi-symmetric  jet.  In  order 
to  be  consistent  with  the  previous  report  the  notations  are  kept 
about  the  same.  An  available  library  plotting  subroutine  is 
utilized  to  present  the  figures  in  the  in-line  plot  so  that  the 
user  can  obtain  the  figures  or  modify  the  scales  without  writing 
additional  subroutines.  Many  figures  are  generated,  but  only  a 
few  are  presented  here  to  demonstrate  the  variations  due  to  the 
changes  of  constants.  In  all  the  figures,  symbols  1,  2,  3 and  4 
represent  the  mean  velocity,  turbulent  shearing  stress,  turbulent 
kinetic  energy  and  dissipation  respectively.  These  quantities 
are  non-dimensionalized  by  the  center  velocity  with  the  different 
scales  indicated.  Figures  3 and  4 show  the  slight  difference  in 
similarity  due  to  two  stations  x-21  and  x*26.  Figure  f>  shows  the 
center  velocity  decay.  It  is  noted  that  the  prediction  is  lower 
than  measurements,  particularly  at  the  development  region. 

However,  further  downstream  the  results  compare  fairly  well  with 
experimenal  evidence.  This  discrepancy  might  be  attributed  to 
the  effects  of  the  initial  conditions,  boundary  layer 
assumptions,  assigned  constants,  or  a small  longitudinal  pressure 
gradient.  It  should  be  enphasized  that  because  of  the  deviation 
in  the  velocity  decay,  the  other  non-dimensional  quantities  such 
as  turbulent  shear,  kinetic  energy,  and  dissipation  are  all 

* BECAUSE  OF  IN-LINE  PLOTTING  * THE  SYMBOLS  4,3»2  AND  1 
MAY  FALL  ON  TOP  OF  EACH  OTHER 
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affected.  It  also  should  be  noted  that  y/x  is  deliberately 
chosen  as  the  coordinate  since  the  y/y^<.  representation  fixes 
one  point  of  y which  make  the  comparison  look  better  than  it 
might  otherwise  have  been.  Figures  6 and  7 show  the  difference 
due  to  one  constant  C$.  Figures  8 and  9 show  the  spread  of 
the  width  due  to  the  entrainment  rate.  From  these  figures,  it 
is  noted  the  results  can  be  different  for  different  constants. 

A comparison  with  experiments  reported  by  Craig  flO]  is  shown 
in  Figure  10.  In  conclusion,  the  results  vary  with  the  values 
of  different  combination  of  constants  and  initial  conditions; 
however,  it  can  compare  fairly  well  with  experiments  provided 
the  constants  are  chosen  properly. 


FURTHER  DEVELOPMENT  AND  SUGGESTIONS 


Since  a computer  program  has  been  satisfactorily  developed  for 
predicting  the  flow  characteristics  such  as  velocity,  turbulent  shear, 
kinetic  energy,  and  dissipation  in  a free  jet,  it  is  recommended  that 
further  investigation  with  more  elaborate  closure  models  to  study 
the  three  components  of  turbulent  intensities  as  well  be  carried  out. 
Initial  investigation  should  focus  on  the  two  dimensional  case 
because  most  of  the  constants  are  determined  by  two  dimensional 
considerations  and  measurements.  After  that,  the  study  should  be 
extended  to  the  axi-symmetric  jet. 


ACKNOWLEDGEMENTS 


Tt»e  author  acknowledges  the  support  from  APOSR  Col.  David  E. 
England,  Col.  Lowell  W.  Ormand  and  Dr.  Brian  Quinn.  The 
author  would  also  like  to  acknowledge  the  opportunity  to  work  at  the 
Ramjet  Technology  Branch,  Air  Force  Aero-Propulsion  Laboratory, 
Wright-Patterson  Air  Force  Base  and  to  express  his  appreciation  for 
the  discussions  with  Dr.  R.R.  Craig  and  Dr.  F.D.  Stull. 


I 

1 


16 


REFERENCES 


1.  Kline,  S.J.,  Morkovin,  M.V.  Sovran,  G. , Cockrell,  D.J. 
"Computation  of  Turbulent  Boundary  Layers"  APOSR-IFP  - 
STANFORD  CONFERENCE  1968,  Department  of  Mechanical  Engineering, 
Stanford  University,  Stanford,  California  94309. 

2.  Reynolds,  W.C. , "Computation  of  Turbulent  Flows  - State-of-the- 
Art,  1970"  Report  MD-27,  Department  of  Mechanical  Engineering, 
Stanford  University,  Stanford,  California  94309. 

3.  Herring,  H.J.  and  Mellor,  G.L.  "A  Method  of  Calculating 
Compressible  Turbulent  Boundary  Layers"  NASA  CR-1444,  1968. 

4.  Patankar,  S.  & Spalding,  D.E. , "Heat  & Mass  Transfer  in  B.L.", 
Inbertext  Books,  1970,  2nd  Edition. 

9.  Hanjalic,  K.  and  Launder,  B.E.,  "A  Reynolds  Stress  Model  of 
Turbulence  and  its  Application  to  Thin  Shear  Flows",  J.  Fluid 
Mech.,  Vol.  92,  Part  4,  pp  609-638,  1972. 

6.  Launder,  B.E.,  "Progress  in  Modeling  of  Turbulent  Transport", 
College  of  Engineering,  Penn  State  University,  Supplementer 
Note  July  28  - August  1,  1975. 

7.  Schlichting,  H. , "Boundary  Layer  Theory”,  6th  Edru  McGraw-Hill, 
N.Y. 

8.  Tsuei,  Y.G. , "Study  of  Reynolds  Stress  Equation  for  Prediction 
of  Flow  Characteristics  of  Free  Jet"  Final  Report  USAF-ASEE 
Summer  Faculty  Research  (WPAFB)  September  1976,  Department 

of  Mechanical  Engineering,  University  of  Cincinnati,  49221. 

9.  Harsha,  P.T. , "Free  Turbulent  Mixing:  A Critical  Evaluation  of 
Theory  and  Experiment".  AEDC-TR-71-36,  ARD  Inc.,  1971. 

10.  Craig,  R.R.,  "Turbulence  In  Free  Diffusion  Flames  of 
Hydrogen-Nitrogen  Mixtures  Burning  In  Still  Air" 
AFAPL-TR-74-68. 


s 

4J 

U 

O 

O 

£ 

r-^ 

(/) 

U-I 

•H 

T5  — 

U 

CNI 

(0 

U-I  VD 

O O' 

J! 

4J  3 

m 

W <D  C , 

C U o 

&)  **H 

w o;  -*-> 

CP  <TJ  1 
>1  C-H 

(U  jC  S 
> 0-0 


w w 

3 *3 

H Cl  U 

«J  w o 

:&S 

£ £3 


^ O (0 
ir . jj  > 
r* 

cr  io  ro 

H H 

o 5 

r—4 

rtj  u-i  ^ r— . 

o • 

4J  aJ  5T 
(D  (U  C 
3 ro  jj 

UH  A)  J 

O)  rtj  W Q 

'Q  > 

c co  ra 
3 a>  c 
(0£  0>w 

J5  u o o 


u ±3 

4J  «*-» 

8 (o 

8S 

c 

<U  4J 

a) 

O'  a> 

sy 

C -H 

§ 

5 c 

> CN 


8 “ 


i5 

.0 

.33 


Initial  Conditions: 
Rend : VI,  U1 , Tl,  ]'J 


boundary  Conditions: 
GU(1)  - GT(1)  - GK(1) 
11U(1)  - 11K(1)  - HE(1) 


CE(1)  - HT(1)  - 0 
1 


Calculate : 

Initial  flow  rate  Q 
Initial  momentum  flux  SM 
Non-dimensional  stream  function  W 


DO  9990  KN  - 1,  KIO: 
This  is  the  main  loo 


X - X + DX 


DO  4 44  ITER  - 1,  NITER 

Iterations,  NITER  - 6 for  KN  £ 3 and  then  NITER  - 3 


DO  120  J * 2,  NP 
Calculation  of  G’s  nnd  H’s 


DO  330  JJ  - 1,  N 

1 

j 

Calculation  of  U,  T,  K,  E 

Divergent 

1 

Stop 

Calculate  entrainment  ENTRN  and 


DO  440  J - 2,  NPP 

Calculate  Y,  Y5.  0 and  momentum  flux  FLUX 


nservation  of  momentum  DM  ■ FLUX  - SM 


anre  DX  as  X Increases 


CALL  OUTPUT 


More 

Run  or  Data? 


o to  the  starting 
point 
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.6.  w r j 


1 


REAL  K f KIfKE 
EXTERNAL  FUNC 

COMMON  LI  ( 25 ) f DU  ( 25 ) f Y < 25 ) f U ( 25  > f E < 25 ) f T < 25 ) rK(  25 ) » TDU  < 25 ) f UDK  ( 25  ) 
COMMON  KNfNMfKKKfNCfNPfNfENTRNfXfQfPSIfSMfFL.UXfDMfYSfRfJRfF 
REAL  UPE  ( 25 ) » GE<25)»  HE(25>»  DY(25)»  YU  (25)  » YUU(25)r  YY(25)f 

1 KK ( 25 ) » DU(25)f  UU ( 25 ) > GU(25)f  HU(25)f  UPU(25)f 

DT  ( 25 ) » TT  ( 25 ) » GT(25)  » HT(25)  r UPT(25)f  UK(25)f  EE(25)f 
U ( 25 ) » riKC25)»  DE  (.25)  t GK(25)»  HK(25)f  UF’K  ( 25 ) » UN  ( 25 ) » 
XX(200)/200>K0.0/r  UC ( 200 ) /200*0 . 0/ f GRAPH ( 400 ) /400*0 . 0/ 
DIMENSION  YI ( 25 ) f UI ( 25 ) f TI ( 25 ) f 01 ( 25 ) f El < 25 ) f KI ( 25 ) 

REAL  YS ( 25 ) f US ( 25 ) f TS ( 25 ) f NS ( 25 ) f ES ( 25 ) 

REAL  XXA<8)/4.f  5.f  A.f  8. f 10. f 15. f 20. f 25./ 

REAL  UCA ( 8 ) /1 . 0 f 0.96f  0.84f  0.70f  0.57f  0.39f  0.29f  0.23/ 

DATA  C/O.O/f  CX  /5.72/fEDGE/10.0/fXG ROW/O. 0025/f  DELTA/0.10/ 

DATA  C1/2.8/fC2/0.45/fCS/0.08/fCEF'S1/1 . 45/ f CEF  S2/2 . 0/ f CEPS/0 . 13/ 
DATA  CC/ 0.07/f  C3/0.00/f  C4/0.00/.  CM/O.OO/f  CONST/O . 8/ f NURITE/1 / 
DATA  FI/3.1 4159/ f PI 1/6 . 28318/ f KNN/10/ f RHO/1 ./ f R JET/1 ./ f UJ/1 . / 
DATA  LL/O/ 

DATA  NPP/15/f YI/O.OOfO. 1 f 0.2f  0.3f  0.4f  0.5f  0.6f  0.7f  .8 f 

1 0 . 85  f 0 . 9 f 0 . 95  f I.Of  I.OSfI.If  10*0.0/ 

DATA  UI/9*1 .Of  0 . 98 f 0.96f0.90f  0.5f  0.25f  O.Of  10*0 . 0/ f JJJ J/2/ 
DATA  TI/O.Of  13*0. OOIf  O.Of  10*0. 0/f  KI/25*0.001/f  EI/25*0.002/ 
DATA  TT ( 1 ) f GU ( 1 ) f GT ( 1 ) f GK ( 1 ) f HT ( l ) f YU ( 1 > fYUU( 1 ) f PRESS/8*0 . 0/ 

DATA  V<1)fUK(1)fYY(1>fGE(1)f  A f DPSI /6*0 . 0/ f XLAST/60 . / f DX/ . 05/ 

DATA  HE < l > f HU ( 1 > fHK< 1 >fUU(1)/4*1.0/f UEDGE/O . / f NNN/10/ fL AMIN/1/ 

CEPS1  = 1 .44 

Cl  = 2.6 

N=NPF— 2 

NP=N+1 

NM=N-1 

NC=NP/2 

C1=C1+0.1 

CS=0.08 

DO  11111  LLLL=1 f 2 
CEPS2=CEPS2-0. 1 
CS=0 . 08 

6666  FORMAT  (/) 

DO  11111  KKKKK=1 f 2 
CS=CS+0.01 

do  mil  jjjjj=1fJjjj 

100  CONTINUE 
ENTRN=0 . 5 
XGROU^-O . 0002 
NNN=1 0 
NWRITE=1 
JR=1 
R=400 • 

RR=R/1000 . 

F=CX* ( C+l . 0 ) 

Y5=l . 


KI  ( NF’P ) =0 . 0 
DO  10  J=lrNPP 
Y(J)=YI(J) 

U( J)=UI( J) 

T(J)=TI(J) 

K ( J ) -KI ( J ) 

E ( J ) =EI < J ) 

10  CONTINUE 

HU ( 1 ) -1 . 0001 
U<1>=0.0 
TDU ( 1 )=0.0 
X = 0.0 
DX-0 . 05 
DM  = 0 4 0 
FLUX=0 . 0 
KKK=68 

WRITE  ( 6 > 6666 ) 

0 = 0.0 
SM=0 . 0 

DO  102  J=2>NPP 
JM= J-l 

DY ( JM  > =Y ( J > -Y ( JM ) 

YRM=0.5*<Y< J)+Y( JM) >*2.*FI 
IF  (JR.EQ.O)  YRM=1 . 

DW ( JM ) = ( U ( J ) +U ( JM ) ) /2 . *DY ( JM ) #YRM 
W( J)=W< JM)+DW< JM> 

0=  0 + 0.5*<U< J)+U( JM) ) *DY(JM)*YRM 
SM=SM+  0.5# (U< J)*#2+  U( JM>**2)*DY< JM)*YRM 

102  CONTINUE 
SMM=SM 

DO  110  J=1»NP 
W( J)=U( J)/W<NPP) 

DW( J)=DW( J)/U(NPP) 

110  CONTINUE 
U(NPP)=1 .0 
PS  1=0 

101  CONTINUE 

DO  9990  KN=1 t KKK 

ENTRN=ENTRN+XGROW*X 

NPPS=NPP 

IF  (KN.LE.3)  NITER=6 
IF  (KN.GT.3)  NITER=3 
XS  = X 
X=X+DX 

DO  103  J=1»NFP 
UK ( J ) =K ( J ) 

103  U(J)=U(J) 

L.AMIN=0 

IF  (LAMIN.EO.l)  R=30 
DO  444  ITER=1» NITER 


DO  120  J=2fNP 
JM~  J- 1 
JF*= J+l 
UEDGE=--0.0 
UCE-LK  1 )-UEDGE 
F=CX*(C  + Y5*(JCE> 

IF  ( (KN.GT.KNN) .AND. (LAMIN.EQ.O) ) F=K(J)/E(J) 

YF=Y( JP)*PII 
YJ=Y< J)*PII 
YM=Y< JM)*PII 
DWPM=U( JP)-W< JM) 

A=0 

B=-ENTRN/F'SI 
RP=RHO/PSI 
FcF'P-4 . 0*RP*RP 

CMU=RPF'/R* ( U ( JM ) *YM*YM+U ( J ) *YJ*YJ ) 

CF’U=RPF'/RJ|<  ( U < JP ) *YP*YF+U  ( J ) *Y  J*Y  J ) 

CMT=CS*RPP*F*<U( JM)*T< JM)*YM*YM+U< J)*T< J)*YJ*YJ> 

CPT=CS*RPP#F*  ( U < JF'  > *T  < JP > * YP*YP+U  ( J ) *T  ( J ) *YJ#YJ  ) 

CMT=CS*RPP*F* ( U ( JM ) *l\ ( JM ) * Yh*  Yh+U ( J ) *K  < J ) *YJ#YJ ) 

CF'T=CS*RPF*F*  ( U ( JP ) *K  ( JP ) *YP* YF  + U ( J ) *l\  ( J ) #Y  J*Y  J ) 

CMK=CONST*CMT 

CPK=CONST*CPT 

CME=0. 5*CEFS*CMT/CS 

CF'E  = 0 . 5*CEPS*CF'T/CS 

G1  = < DU  < J > /DX+4 . *A+B* < U ( JP ) +3 . *U ( J ) ) ) /A . /DUPM 
G2=(3»/DX-B)/4. 

G3=  ( DU ( JM ) /DX-4 . *A-B* ( U ( JM ) +3 . *U ( J ) ) ) /A . /DUPM 
G5U=CPU/DWPM/DW ( J ) 

G5T=CPT /DUPM/DU < J ) 

G5K=CPK/DWPM/DW < J ) 

G5E=CPE/DUPM/DU ( J ) 

G6U=CMU/DUPM/DU ( JM ) 

G6T=CMT /DUF’M/DU  ( JM ) 

G6K=CMK/DUPM/DU ( JM ) 

G6E=CME/DUPM/DU ( JM ) 

AU=G3-G6U 

AT=G3-G6T 

AK=G3-G6K 

AE=G3 -G<6E 

BU=G2+G5(J+G6U 

BT =G2+G5T +G6T  + C1/F/U ( J ) 

BK=G2+G5K+G6K+1 . /F/U ( J ) 

BE=G2+G5E+G6E+CEPS2/F/U< J> 

CU=G1-G5U 
CT=G1-G5T 
CK=G1 -G5K 
CE=G1 -G5E 

IF  (ITER.GT.l)  GO  TO  111 

UF'U  < J ) = < DW  ( JM ) *U  ( JM ) +3 . *DUPM*U  < J ) +DU  ( J ) *U  ( JP ) ) /A . /DX/DUPM 
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UF'T  < J ) = ( DU  ( JM ) *T  ( JM  ) +3  . *DWPM*T  ( J ) + DU  ( J ) *T  ( JP ) ) /A . /DX/DUPM 
UPK ( J ) = ( DU ( JM ) *K ( JM ) +3 . *DUPM*K ( J ) + DU < J > *K ( JP  > ) /A . /DX/DUPM 
LIFE ( J ) = ( DU ( JM ) *E  < JM ) +3 . *DUPM*E ( J ) +DU ( J ) *E ( JP ) ) /A . /DX/DUPM 
111  CONTINUE 

DUDU= ( U ( JP ) -U  < JM ) ) /DUPM 
DTDU= ( T ( JP ) *YF-T ( JM ) * YM ) /DUPM 
T JJM=  ( T ( J ) -T  < JM  ) ) /DU  < JM ) * ( Y J+YM  ) /2 , 

T JP  J= ( T ( JP ) -T  < J ) ) /DU  < J ) * ( YF'  + Y J ) /2 . 

DTDU=  < T JP J+T J JM ) /2 . 

SU=FRESS-RP*DTDU 

ST  = -C2*RP*K< J)*DUDW  *YJ 

SK=  -RP*T< J)*DUDW*YJ 

SE=-CEPS1*RF’*T(  J)/F#DUDU*YJ 

FU=UPU ( J ) +SU 

FT =UF'T  ( J ) +ST 

FK=UPK ( J ) +SK 

FE  = UF'E  ( J ) +SE 

GU(J)=  (FU-AU*GU< JM) ) / ( BU+AU*HU ( JM ) ) 

HU( J)=-CU/<BU+AU*HU< JM) ) 

GT( J)=(FT-AT*GT( JM) > / < BT+AT*HT < JM) ) 

FIT ( J ) =-CT  /(BT+AT*HT<  JM)  ) 

GK< J)=(FK-AK*GK( JM) >/<BK+AK*HK< JM) ) 

HK(J)=-CK  /(BK+AK*HK'(  JM)  ) 

GE(  J)  = (FE-AE*GE( JM) ) / ( BE  + AE*HE ( JM ) ) 

HE( J)=-CE/(BE+AE*HE( JM) ) 

120  CONTINUE 


130 

230 


330 


435 


DO  330  JJ=lfNP 

U ( NPF'-J  J ) =HU  ( NPP-JJ ) *U  ( NPP- J J+ 1 ) +GU  ( NPP- J J ) 

T ( NPP-JJ ) =HT ( NPP-JJ ) *T ( NPP- J J+ 1 ) +GT ( NPP-JJ ) 

K ( NPP-JJ ) =HK  C NPP-JJ ) *K ( NPP- J J+ 1 ) +GK < NPP-JJ ) 

E (NPP-JJ ) =HE ( NPP-JJ ) *E ( NPP- J J+l ) +GE ( NPP-JJ ) 

CONTINUE 

U < NPP )=U(NP> /EDGE 
T < NPP >=T<NP> /EDGE 
K < NPP >=K(NP> /EDGE 
E ( NPP )=E(NF> /EDGE 
0=0.0 
FLUX=0 . 0 
DO  440  J=2»NPP 
JM= J-l 

YRM=0.5*(Y( J)+Y( JM) )*2.*PI 
IF  (JR.EQ.O)  YRM=1 ♦ 

DY ( JM ) =2 . *PSI *DU ( JM ) / ( U ( J ) +U  < JM ) ) /YRM 
Y ( J ) =Y ( JM ) +DY ( JM ) 

0=  0 + 0.5*(U( J)+U< JM) ) *DY(JM)*YRM 
FLUX=FLUX  + 0.5*(U( J)**2+  U ( JM ) **2 > *DY ( JM ) *YRM 
UU ( J ) = ( U ( J ) -UEDGE ) /UCE 

IF  < . NOT .(UU(J).LT.O.S. AND .UU(JM).GT.0.5))  GO  TO  435 
Y5=Y( JM)+(0.5-UU( JM) )/(UU< J)-UU( JM) >*(Y( J)-Y( JM) ) 
CONTINUE 


IF  (J.EG.NPP)  GO  TO  440 
JF-J+1 

DWPM=W( JF)-W( JM) 

ENTU=ENTRN  *W ( J ) /PSI * ( U ( JP ) -U ( JM ) ) /DWPM 
ENTK=ENTRN  *W ( J ) /PSI* < K ( JP ) -K ( JM ) ) /DWPM 
TDU (J)=U(J)#((U(J)-V(J)) /DX-ENTU  > 

UDK ( J ) =U ( J ) * ( ( K < J ) -VK C J ) ) /DX-ENTK ) * Y5/UCE**3 
440  CONTINUE 
DM=FLUX-SM 

IF  ( (UCNP) .GT.U(N) > .AND. (T(NP) .GT.T(N)  ) ) GO  TO  9999 
444  CONTINUE 

UC(KN)=U< 1 > 

XX ( KN ) = X/2 . 

IF  (XX(KN) .LE.1.1 > XX ( KN ) = 1 . 1 

DPSI=ENTRN*DX 

PSI=F'SI+DPSI 

IF  (KN.EG.l)  WRITE  (6*1100) 

1100  FORMAT  (IX*'  KN  X PSI  YNPP  Y5  ENTRN  DM  U(l>  U<3>  U 

1(4)  U ( 5 ) U ( 6 ) UNM  UN  UNP  T(NC>  T ( N ) K(l)  KC  KNM 

2KN  KNP ' ) 

IF  ( (KN.LE.l) .OR. (KN.EQ.NWRITE)  ) CALL  OUTPUT  (5) 

IF  (KN.GT.51)  NNN=5 

IF  (KN.EQ.NWRITE)  NWRITE=NWRITE+NNN 
480  CONTINUE 
490  CONTINUE 

IF  (X.GE.1.0)  DX=DELT A*X 
IF  (DX.GT.2.0)  DX=2 . 0 

IF  ( JJJJJ.EQ.2. AND.X.LE.5,0)  GO  TO  556 

DO  555  J=1*NPP 

U( J)=U( J)*SQRT(SM/FLUX) 

555  CONTINUE 

556  CONTINUE 
SMM=FLUX 

IF  ( .NOT. ( (KN.EQ.  63 ) . OR . ( KN . EQ . 82 ) . OR . ( KN . EQ . KKK ) > ) GO  TO  9990 
DO  1 J=1»NP 
GRAPH ( J)=Y( J)/X 

IF  ( GRAPH (J).GT.0.16>  GRAPH ( J ) =0 . 16 
GRAPH ( J+NP )=U(J)/U(1) 

GRAPH ( J+2#NP )=T(J)/U(1) ##2#40 . 

GRAPH ( J+3*NP )=K( J)/U(l ) **2*10 . 

GRAPH  ( .J+4*NP  )=E(  J)  *Y5/U  ( 1 ) **3*40 . 

1 CONTINUE 

CALL  PL0T9  ( 1 * GRAPH  * NF'  * 5 * 0 * FUNC  * .2  * 0.0*1. 4*  0.*  400) 

CALL  OUTPUT  (6) 

WRITE  (6*7770) 

WRITE  (6*7771)  JR*  N*  KKK  * EDGE  * CEPS1,  Q *SM*FLUX*  DM*  X * F*  C* 
1 CX»  Cl*  C2*  CS*  CONST*  DX  *RR* ENTRN 
WRITE  (6*7772)  CEPS*CEPS2 
WRITE  (6*6666) 

CALL  OUTPUT  (5) 


I 
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9990  CONTINUE 
KNM=KN-1 
KNM8=KNM+8 
DO  2 J=lrKNM 
GRAPH <J+8)=XX(J) 

2 GRAPH ( J+KNM8+8 ) =UC ( J) 

DO  3 J=1 » 8 

GRAPH ( J ) =XXA ( J ) 

3 GRAPH  < J+KNM8 ) =UCA ( J ) 

CALL  PL0T9  ( 1 , GRAPH , KNM8 » 2 » 0 » FUNC  , 28. r 0.0»1.4r  O.f  400) 

KNP=I\N  + 1 
9999  CONTINUE 
700  WRITE  (.6,7770) 

7770  FORMAT  (/,'  JR  N NKK  EUGE  CEPS1  0 SM  FLUX  HM 

1 X F C CX  Cl  C2  CS  CONST  HX  R/ 

21000  ENTRN  ' ) 

RR=R/1000. 

WRITE  (6,7771)  JR,  M,  KKK , EDGE , CEP SI*  Q f SM r FLUX , Dh,  X , F,  C, 
1 CX»  Cl,  C2r  CS,  CONST , DX  , RR , ENTRN 

7771  FORMAT  ( IX » 31 4 , 17F7 . 2 ) 

WRITE  (6,6666) 

WRITE  (6,7772)  CEPSrCEPS2 

7772  FORMAT  ( IX r ' CEPS= ' , F5 . 3 » ' CEPS2= ' , F5 . 3 > 

11111  CONTINUE 

STOP 

END 
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SUBROUTINE  OUTPUT  (ICALL) 

REAL  K»  KKt  KI 

COMMON  W < 25 ) f DW < 25  > f Y < 25 > f U < 25  > f E < 25  > f T ( 25  > f K < 25 ) f THU (25) f UDK ( 25 ) 
COMMON  KN  f NM  f KKK  f NC  f NP  f N f ENTRN  f X f Q f PS I f SM  f FLUX  f DM  f Y5  f R f JR  f F 
DIMENSION  YY(25)fTT(25)fKK(25)fEE(25)fUU(25) 

NPF-N+2 

1111  FORMAT  <I4f1P11E10.2f  1P2E9.1) 

UELGE=0 » 0 

UCE=U< 1 )-UEDGE 
UCEN=UCE/ ( 1 . O-UEDGE ) 

Z=X/2 . 

GO  TO  <100f200f200f200f500f600) f ICALL 
100  WRITE  (6f 1000) 

1000  FOf  1AT  <1Xf'  INITIAL  AND  BOUNDARY  CONDITIONS  FOR  NPP  W DW  Y 
1 U T K ' > 

1112  FORMAT  < IX  f 14  f 25F5 . 2 ) 

WRITE  <6f1112>  NPPf  < W < J > f J=1 f NPP > 

WRITE  < 6 f 1 1 12 ) N f (DW( J> f J=1 fNF) 

WRITE  (6f1112)  NPPf  ( Y ( J ) f J=1 f NPP ) 

WRITE  <6f1112>  NPPf  < U ( J > f J=1 f NPP ) 

WRITE  (6f1112>  NPFf  (T< J) f J=1 fNPF) 

WRITE  <6f 1112)  NPFf  (K( J) f J=1 fNPP) 

WRITE  ( 6 f 6666 > 

6666  FORMAT  </> 

200  RETURN 

500  WRITE  ( 6 f 1234 ) KNfZfFSIf  Y(NPP)f  Y5f  ENTRNf  DMf  UCENfU(3)f  U(4)f 

1 U(5) f U ( 6 ) f U(NM) f U(N)f  U(NP) fT(NC) fT(N) f K(1)f  K(NC)f 

2 K ( NM ) f K ( N ) f N(NF') 

1234  FORMAT  ( I4f 14F6.2f7F6.3> 

RETURN 

600  WRITE  ( 6 f 6660 ) 

6660  FORMAT  </f1Xf'  J W Y U T K 

1 RX  YY  UU  TT  KK  EE  TDU 

2 UDK  ' ) 

DO  610  J=2  f NP  f 2 
RX=Y( J>/X 
UU( J>=U( J>/UCE 
YY( J)=Y ( J)/Y5 
TT< J>=T< J)/UCE**2 
KK( J)=K( J)/UCE**2 
EE< J>=E( J)*Y5/UCE**3 

610  WRITE  (6f1111)  JfW(J)f  Y(J)f  U(J)f  T<J)fK(J)f  RX  f YY(J)f  UU<J>f 
1 TT ( J) f KK ( J) f EE < J ) f TDU(J)f  UDK(J) 

WRITE  ( 6 f 6666 ) 

RETURN 

END 


SUBROUTINE  FUNC  (XINfYOUT) 
Y0UT=0. 

RETURN 

END 
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FIG. 3 PROFILES  OF  VELOCITY.  SHEAR . ENERGY  AND  DISSIPATION 
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FIG. 7 PROFILES  OF  VELOCITY,  SHEAR,  ENERGY  AND  DISSIPATION 
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FIG. 8 PROFILES  OF  VELOCITY*  SHEAR  * ENERGY  AND  DISSIPATION 
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COMPARISON  WITH  EXPERIMENTS 
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